

# 07-05-2021
# REPLICATION FILES

# Commodity shocks and incumbency effects - Novaes & Schiumerini

#This file replicates the following figures

# Main text figures : 1, 2, 3
# Appendix: incumbency effects with controls, scatter plot figure incumbency Effects; density; density tests, histogram
# You should Restart R before anything

# These scripts have been run in the following system and R version:

# platform       x86_64-apple-darwin17.0     
# arch           x86_64                      
# os             darwin17.0                  
# system         x86_64, darwin17.0          
# status                                     
# major          4                           
# minor          0.2                         
# year           2020                        
# month          06                          
# day            22                          
# svn rev        78730                       
# language       R                           
# version.string R version 4.0.2 (2020-06-22)
# nickname       Taking Off Again   


required_packages = c('estimatr','data.table','rdd','apsrtable','stargazer','rdrobust','multiwayvcov','lmtest','ggplot2','ggthemes','rddensity')
lapply(required_packages, require, character.only = TRUE)

library(data.table)
library(ggplot2)
library(ggthemes)
library(gridExtra)
library(sandwich)
library(rdd)
library(rddensity)

setwd(dir = '.') # Please insert folder where files are stored

load("replication_data.RData")


# Figure 1: Incumbency effects plot in the appendix ----

rd <- subset(retro, (year <2016 & year>2000) & ineligible2 == 0)
# ineligible2 == 0 removes municipalities with incumbents that would be term-limited at t+1 if they win at t

rd[,y2000:=ifelse(year==2000,1,0)]
rd[,y2004:=ifelse(year==2004,1,0)]
rd[,y2008:=ifelse(year==2008,1,0)]
rd[,y2012:=ifelse(year==2012,1,0)]

all = rdrobust(y=rd$f_elected, x = rd$margin, all = TRUE, kernel = 'triangular', covs = cbind(rd$y2012, rd$y2008), cluster = cbind(rd$muni_code))

rd_rural <- rd[ruralPop>median(ruralPop,na.rm=TRUE)]
rural = with(rd_rural, rdrobust(y = f_elected, x = margin, all=TRUE, kernel = 'triangular', covs = cbind(y2012, y2008),cluster = cbind(muni_code)))

rd_urban <- rd[ruralPop < median(ruralPop,na.rm=TRUE)]
urban = with(rd_urban, rdrobust(y = f_elected, x = margin, all=TRUE, kernel = 'triangular', covs = cbind(y2012, y2008), cluster = cbind(muni_code)))


late.plot = data.frame(group = c('All','Rural','Urban'),
                       late = c(all$coef[[3]],rural$coef[[3]],
                                urban$coef[[3]]),
                       se = c(all$se[[3]],
                              rural$se[[3]],urban$se[[3]]),
                       ci_u = c(all$ci[[3,1]],
                                rural$ci[[3,1]],
                                urban$ci[[3,1]]),
                       ci_l = c(all$ci[[3,2]],
                                rural$ci[[3,2]],
                                urban$ci[[3,2]]))

incumbency = ggplot(late.plot,aes(group=group,x=group)) +
  geom_point(aes(y=late*100,x=group),size=2) +
  geom_hline(yintercept = 0,linetype = 'dashed') +
  theme_minimal() +
  ylab('Probability Win Next Election (in p.p.)') +
  xlab('Group of Municipalities') +
  geom_errorbar(aes(ymax = ci_u*100, ymin = ci_l*100), width=0,size=.5) 

#ggsave('incumbency.pdf', plot = incumbency, device = 'pdf',height = 10, width = 10, units = 'cm')

# Figure 2, rural municipalities ----

rd <- subset(retro, (year < 2016 & year > 2000) & ineligible2 == 0)

rd = rd[ruralPop > median(retro$ruralPop,na.rm=TRUE)] # Rural 
bw.u <- rdbwselect(rd$f_elected,rd$margin)
bwselect = bw.u$bws[[2]] * 100 # The line in the plot

rdA = rd[inflation >= (median(rd$inflation,na.rm=T))] # Above median inflation
rdB = rd[inflation < (median(rd$inflation,na.rm=T))]  # Below median inflation

#plot parameters
xlim = c(0,15)
ylim= c(-.6,.3)
xlimS = c(-.15,.15)
ylimS= c(0,.5)

#Sequence
sequence=seq(0.01, .15, .01)

#Above median
tt=rdA
tt[,dv:= f_elected]

data.plot.1 <- data.frame(coef = as.numeric(),
                          SD = as.numeric(),
                          n = as.numeric(),
                          stringsAsFactors = FALSE)

for(i in sequence){
  rd.data.plot.1 <- tt[abs(tt$margin) <= i, ]
  weights <- kernelwts(rd.data.plot.1$margin, center = 0, bw = i, kernel= 'triangular')
  data.plot.1reg <- with(rd.data.plot.1, lm(I(dv) ~ elected*margin + as.factor(year)))
  data.plot.1reg.SE <- sqrt(diag(cluster.vcov(data.plot.1reg, cbind(rd.data.plot.1$muni_code))))
  tempo.a <- data.frame(coef = data.plot.1reg$coefficients[2],
                        SD = data.plot.1reg.SE[2],
                        n = data.plot.1reg$df)
  data.plot.1 <- rbind(data.plot.1,tempo.a)   
}

data.plot.1$cuts <- sequence
data.plot.1$confidence <- with(data.plot.1, qt(.975, df = n - 1) * SD)
data.plot.1$upper <- with(data.plot.1, coef + confidence)
data.plot.1$lower <- with(data.plot.1, coef - confidence)

data.plot.1$description <- 'Probability Winning - Rural Municipalities'
data.plot <- rbind(data.plot.1)
data.plot_above <- copy(data.plot)

plot.above <- ggplot(data.plot, aes(x = cuts * 100, y = coef, group = description)) +
  geom_line(aes(color = description), size = .5, alpha = .7) +
  geom_line(mapping = aes(y = upper, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_line(mapping = aes(y = lower, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_point(data = data.plot, mapping= aes(x=cuts * 100, y=coef),alpha=.9 ,na.rm=T,size=2) +
  geom_hline(yintercept=0, linetype="dotted",colour='gray10') +
  geom_vline(xintercept = bwselect, colour='gray10', alpha=.5) +
  theme_minimal() +
  xlab("") +
  ylab('(1.B)\nLATE P(win, t+1)') +
  scale_color_grey(name="") +
  annotate('text', y = -0.5, x = 9.5, label = 'optimal\nbandwidth',size=4,col='gray40',fontface = 'italic',family='Times') +
  coord_cartesian(ylim=ylim,xlim=xlim) +
  theme(legend.position="none") +
  geom_text(aes(y = 0.1,label = n),vjust=-.50,size=2,hjust=-.1,col='gray40',angle=60)

dos <- subset(tt, margin > 0)
uno <- subset(tt, margin < 0)
tres <- rbind(uno,dos)

y = (tres$dv)

cut <- cut(tres$margin,400, include.lowest = TRUE)
tmp <- aggregate(y, by=list(cut = cut), FUN=mean, na.rm=T)
tmp1 <- aggregate((tres$margin), by=list(cut = cut), FUN=mean, na.rm=T)
data <- data.frame(margin = tmp1$x, y = tmp$x)

scatterAbove = ggplot() +
  stat_smooth(data = tres, aes(margin,y,group=elected), method="loess",na.rm=T,se=FALSE, level = .95, fill = 'gray80',colour="black", size=1.25 , formula = y ~ x) +
  geom_point(data = data[data$margin <0,], aes(margin, y), alpha=.9 ,na.rm=T,size=2) +
  geom_point(data = data[data$margin >0,], aes(margin, y), alpha=.9 ,na.rm=T,size=2,shape = 18) +
  theme_minimal() +
  theme(plot.margin=unit(c(1,0,0,1),"mm")) +
  xlab("") +
  geom_vline(xintercept=0, linetype="dotted") +
  coord_cartesian(ylim=ylimS,xlim=xlimS) +
  scale_y_continuous(labels=scales::percent) +
  scale_x_continuous(labels=scales::percent) +
  ylab("(1.A)\nProbability Win Next Election")

allA = arrangeGrob(scatterAbove,plot.above,nrow=1,top = '(1) Positive Shock\n(above median)')


tt=rdB
tt[,dv:= f_elected]

data.plot.1 <- data.frame(coef = as.numeric(),
                          SD = as.numeric(),
                          n = as.numeric(),
                          stringsAsFactors = FALSE)
for(i in sequence){
  rd.data.plot.1 <- tt[abs(tt$margin) <= i, ]
  weights <- kernelwts(rd.data.plot.1$margin, center = 0, bw = i, kernel= 'triangular')
  data.plot.1reg <- with(rd.data.plot.1, lm(I(dv) ~ elected*margin+as.factor(year), weights = weights))
  data.plot.1reg.SE <- sqrt(diag(cluster.vcov(data.plot.1reg, cbind(rd.data.plot.1$muni_code))))
  tempo.a <- data.frame(coef = data.plot.1reg$coefficients[2], SD = data.plot.1reg.SE[2], n = data.plot.1reg$df)
  data.plot.1 <- rbind(data.plot.1,tempo.a)   
}


data.plot.1$cuts <- sequence
data.plot.1$confidence <- with(data.plot.1, qt(.975, df = n - 1) * SD)
data.plot.1$upper <- with(data.plot.1, coef + confidence)
data.plot.1$lower <- with(data.plot.1, coef - confidence)


data.plot.1$description <- 'Probability Winning - Rural Municipalities'
data.plot <- rbind(data.plot.1)
data.plot_below <- copy(data.plot)

plot.below <- ggplot(data.plot, aes(x = cuts * 100, y = coef, group = description)) +
  geom_line(aes(color = description), size = .5, alpha = .7) +
  geom_line(mapping = aes(y = upper, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_line(mapping = aes(y = lower, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_point(data = data.plot, mapping= aes(x=cuts * 100, y=coef),alpha=.9 ,na.rm=T,size=2) +
  geom_hline(yintercept=0, linetype="dotted",colour='gray10') +
  geom_vline(xintercept=bwselect, colour='gray10',alpha=.5) +
  theme_minimal() +
  xlab("Absolute Distance\nTreatment and Control (p.p.)") +
  ylab('(2.B)\nLATE P(win, t+1)') +
  scale_color_grey(name="") +
  annotate('text', y = -0.5, x = 9.5, label = 'optimal\nbandwidth',size=4,col='gray40',fontface = 'italic',family='Times') +
  coord_cartesian(ylim=ylim,xlim=xlim) +
  theme(legend.position="none") +
  geom_text(aes(y = 0.1,label = n),vjust=-.50,size=2,hjust=-.1,col='gray40',angle=60)


dos <- subset(tt, margin > 0)
uno <- subset(tt, margin < 0)
tres <- rbind(uno,dos)

y = (tres$dv)

cut <- cut(tres$margin,400, include.lowest = TRUE)
tmp <- aggregate(y, by=list(cut = cut), FUN=mean, na.rm=T)
tmp1 <- aggregate((tres$margin), by=list(cut = cut), FUN=mean, na.rm=T)
data <- data.frame(margin = tmp1$x, y = tmp$x)

scatterBelow = ggplot() +
  stat_smooth(data = tres, aes(margin,y,group=elected), method="loess", na.rm=T, se=FALSE, level = .95, fill = 'gray80',colour="black", size=1.25, formula = y ~ x) +
  geom_point(data = data[data$margin <0,], aes(margin, y),alpha=.9 ,na.rm=T,size=2) +
  geom_point(data = data[data$margin >0,], aes(margin, y),alpha=.9 ,na.rm=T,size=2,shape = 18) +
  theme_minimal() +
  theme(plot.margin=unit(c(1,0,0,1),"mm")) +
  xlab("Margin") +
  geom_vline(xintercept=0, linetype="dotted") +
  coord_cartesian(ylim=ylimS,xlim=xlimS) +
  scale_y_continuous(labels=scales::percent) +
  scale_x_continuous(labels=scales::percent) +
  ylab("(2.A)\nProbability Win Next Election")

allB = arrangeGrob(scatterBelow,plot.below,nrow=1,top = '(2) Negative Shock\n(below median)')

plot = grid.arrange(allA,allB)

#ggsave('scatterplot_rural.pdf', plot = plot, device = 'pdf',width = 24, height = 20, units = 'cm')

# Figure 3, urban municipalities ----

rd <- subset(retro, (year < 2016 & year > 2000) & ineligible2 == 0)

rd = rd[ruralPop < median(retro$ruralPop,na.rm=TRUE)] # Urban 
bw.u <- rdbwselect(rd$f_elected,rd$margin)
bwselect = bw.u$bws[[2]] * 100 # The line in the plot

rdA = rd[inflation >= (median(rd$inflation,na.rm=T))] # Above median inflation
rdB = rd[inflation < (median(rd$inflation,na.rm=T))]  # Below median inflation

#plot parameters
xlim = c(0,15)
ylim= c(-.6,.3)
xlimS = c(-.15,.15)
ylimS= c(0,.5)

#Sequence
sequence=seq(0.01, .15, .01)

#Above median
tt=rdA
tt[,dv:= f_elected]

data.plot.1 <- data.frame(coef = as.numeric(),
                          SD = as.numeric(),
                          n = as.numeric(),
                          stringsAsFactors = FALSE)

for(i in sequence){
  rd.data.plot.1 <- tt[abs(tt$margin) <= i, ]
  weights <- kernelwts(rd.data.plot.1$margin, center = 0, bw = i, kernel= 'triangular')
  data.plot.1reg <- with(rd.data.plot.1, lm(I(dv) ~ elected * margin + as.factor(year)))
  data.plot.1reg.SE <- sqrt(diag(cluster.vcov(data.plot.1reg, cbind(rd.data.plot.1$muni_code))))
  tempo.a <- data.frame(coef = data.plot.1reg$coefficients[2],
                        SD = data.plot.1reg.SE[2],
                        n = data.plot.1reg$df)
  data.plot.1 <- rbind(data.plot.1,tempo.a)   
}

data.plot.1$cuts <- sequence
data.plot.1$confidence <- with(data.plot.1, qt(.975, df = n - 1) * SD)
data.plot.1$upper <- with(data.plot.1, coef + confidence)
data.plot.1$lower <- with(data.plot.1, coef - confidence)

data.plot.1$description <- 'Probability Winning - Urban Municipalities'
data.plot <- rbind(data.plot.1)
data.plot_above <- copy(data.plot)

plot.above <- ggplot(data.plot, aes(x = cuts * 100, y = coef, group = description)) +
  geom_line(aes(color = description), size = .5, alpha = .7) +
  geom_line(mapping = aes(y = upper, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_line(mapping = aes(y = lower, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_point(data = data.plot, mapping= aes(x=cuts * 100, y=coef),alpha=.9 ,na.rm=T,size=2) +
  geom_hline(yintercept=0, linetype="dotted",colour='gray10') +
  geom_vline(xintercept = bwselect, colour='gray10', alpha=.5) +
  theme_minimal() +
  xlab("") +
  annotate('text', y =-0.4, x = 13, label = 'optimal\nbandwidth',size=4,col='gray40',fontface = 'italic',family='Times') +
  ylab('(1.B)\nLATE P(win, t+1)') +
  scale_color_grey(name="") +
  coord_cartesian(ylim=ylim,xlim=xlim) +
  theme(legend.position="none") +
  geom_text(aes(y = 0.1,label = n),vjust=-.50,size=2,hjust=-.1,col='gray40',angle=60)


dos <- subset(tt, margin > 0)
uno <- subset(tt, margin < 0)
tres <- rbind(uno,dos)

y = (tres$dv)

cut <- cut(tres$margin,400, include.lowest = TRUE)
tmp <- aggregate(y, by=list(cut = cut), FUN=mean, na.rm=T)
tmp1 <- aggregate((tres$margin), by=list(cut = cut), FUN=mean, na.rm=T)
data <- data.frame(margin = tmp1$x, y = tmp$x)

scatterAbove = ggplot() +
  stat_smooth(data = tres, aes(margin,y,group=elected), method="loess",na.rm=T,se=FALSE, level = .95, fill = 'gray80',colour="black", size=1.25 , formula = y ~ x) +
  geom_point(data = data[data$margin <0,], aes(margin, y), alpha=.9 ,na.rm=T,size=2) +
  geom_point(data = data[data$margin >0,], aes(margin, y), alpha=.9 ,na.rm=T,size=2,shape = 18) +
  theme_minimal() +
  theme(plot.margin=unit(c(1,0,0,1),"mm")) +
  xlab("") +
  geom_vline(xintercept=0, linetype="dotted") +
  coord_cartesian(ylim=ylimS,xlim=xlimS) +
  scale_y_continuous(labels=scales::percent) +
  scale_x_continuous(labels=scales::percent) +
  ylab("(1.A)\nProbability Win Next Election")

allA = arrangeGrob(scatterAbove,plot.above,nrow=1,top = '(1) Positive Shock\n(above median)')

# Negative

tt=rdB
tt[,dv:= f_elected]

data.plot.1 <- data.frame(coef = as.numeric(),
                          SD = as.numeric(),
                          n = as.numeric(),
                          stringsAsFactors = FALSE)
for(i in sequence){
  rd.data.plot.1 <- tt[abs(tt$margin) <= i, ]
  weights <- kernelwts(rd.data.plot.1$margin, center = 0, bw = i, kernel= 'triangular')
  data.plot.1reg <- with(rd.data.plot.1, lm(I(dv) ~ elected*margin+as.factor(year), weights = weights))
  data.plot.1reg.SE <- sqrt(diag(cluster.vcov(data.plot.1reg, cbind(rd.data.plot.1$muni_code))))
  tempo.a <- data.frame(coef = data.plot.1reg$coefficients[2],
                        SD = data.plot.1reg.SE[2],
                        n = data.plot.1reg$df)
  data.plot.1 <- rbind(data.plot.1,tempo.a)   
}


data.plot.1$cuts <- sequence
data.plot.1$confidence <- with(data.plot.1, qt(.975, df = n - 1) * SD)
data.plot.1$upper <- with(data.plot.1, coef + confidence)
data.plot.1$lower <- with(data.plot.1, coef - confidence)


data.plot.1$description <- 'Probability Winning - Urban Municipalities'
data.plot <- rbind(data.plot.1)
data.plot_below <- copy(data.plot)

plot.below <- ggplot(data.plot, aes(x = cuts * 100, y = coef, group = description)) +
  geom_line(aes(color = description), size = .5, alpha = .7) +
  geom_line(mapping = aes(y = upper, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_line(mapping = aes(y = lower, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_point(data = data.plot, mapping= aes(x=cuts * 100, y=coef),alpha=.9 ,na.rm=T,size=2) +
  geom_hline(yintercept=0, linetype="dotted",colour='gray10') +
  geom_vline(xintercept=bwselect, colour='gray10',alpha=.5) +
  theme_minimal() +
  xlab("Absolute Distance\nTreatment and Control (p.p.)") +
  ylab('(2.B)\nLATE P(win, t+1)') +
  scale_color_grey(name="") +
  coord_cartesian(ylim=ylim,xlim=xlim) +
  theme(legend.position="none") +
  annotate('text', y =-0.4, x = 13, label = 'optimal\nbandwidth',size=4,col='gray40',fontface = 'italic',family='Times') +
  geom_text(aes(y = 0.1,label = n),vjust=-.50,size=2,hjust=-.1,col='gray40',angle=60)


dos <- subset(tt, margin > 0)
uno <- subset(tt, margin < 0)
tres <- rbind(uno,dos)

y = (tres$dv)

cut <- cut(tres$margin,400, include.lowest = TRUE)
tmp <- aggregate(y, by=list(cut = cut), FUN=mean, na.rm=T)
tmp1 <- aggregate((tres$margin), by=list(cut = cut), FUN=mean, na.rm=T)
data <- data.frame(margin = tmp1$x, y = tmp$x)

scatterBelow = ggplot() +
  stat_smooth(data = tres, aes(margin,y,group=elected), method="loess", na.rm=T, se=FALSE, level = .95, fill = 'gray80',colour="black", size=1.25, formula = y ~ x) +
  geom_point(data = data[data$margin <0,], aes(margin, y),alpha=.9 ,na.rm=T,size=2) +
  geom_point(data = data[data$margin >0,], aes(margin, y),alpha=.9 ,na.rm=T,size=2,shape = 18) +
  theme_minimal() +
  theme(plot.margin=unit(c(1,0,0,1),"mm")) +
  xlab("Margin") +
  geom_vline(xintercept=0, linetype="dotted") +
  coord_cartesian(ylim=ylimS,xlim=xlimS) +
  scale_y_continuous(labels=scales::percent) +
  scale_x_continuous(labels=scales::percent) +
  ylab("(2.A)\nProbability Win Next Election")

allB = arrangeGrob(scatterBelow,plot.below,nrow=1,top = '(2) Negative Shock\n(below median)')

plot = grid.arrange(allA,allB)

#ggsave('scatterplot_urban.pdf', plot = plot, device = 'pdf',width = 24, height = 20, units = 'cm')

#Figure in appendix, IE with controls

retro[, gdp_pc := gdp / pop]
retro[, pbf_pc := total_pbf_families / pop]
retro[, avg_age := mean(age), by = .(muni_code, year)]
retro[, avg_experience := mean(experience), by = .(muni_code, year)]
retro[, total_left := sum(leftwing), by = .(muni_code, year)]
retro[, total_cands := .N, by = .(muni_code, year)]

illit <- unique(retro[!is.na(illiteracy), .(muni_code, year, illiteracy)])
setkey(illit, muni_code)
setkey(retro, muni_code)
retro$illiteracy = NULL
retro[illit, illiteracy := illiteracy]

rd <- subset(retro, (year <2016 & year>2000) & ineligible2 == 0)

rd[,y2000:=ifelse(year==2000,1,0)]
rd[,y2004:=ifelse(year==2004,1,0)]
rd[,y2008:=ifelse(year==2008,1,0)]
rd[,y2012:=ifelse(year==2012,1,0)]

all = rdrobust(y=rd$f_elected, x = rd$margin, all = TRUE, kernel = 'triangular', 
               covs = cbind(rd$y2012, rd$y2008,rd$gdp_pc, rd$nonwhite, rd$gini, rd$avg_age, rd$illiteracy, rd$avg_experience, rd$total_left, rd$total_cands, rd$pbf_pc), 
               cluster = cbind(rd$muni_code))

rd_rural <- rd[ruralPop>median(ruralPop,na.rm=TRUE)]
rural = with(rd_rural, rdrobust(y = f_elected, x = margin, all=TRUE, kernel = 'triangular', covs = cbind(y2012, y2008,gdp_pc, nonwhite, gini, avg_age, illiteracy, avg_experience, total_left, total_cands, pbf_pc),cluster = cbind(muni_code)))

rd_urban <- rd[ruralPop < median(ruralPop,na.rm=TRUE)]
urban = with(rd_urban, rdrobust(y = f_elected, x = margin, all=TRUE, kernel = 'triangular', covs = cbind(y2012, y2008,gdp_pc, nonwhite, gini, avg_age, illiteracy, avg_experience, total_left, total_cands, pbf_pc), cluster = cbind(muni_code)))


late.plot = data.frame(group = c('All','Rural','Urban'),
                       late = c(all$coef[[3]],rural$coef[[3]],
                                urban$coef[[3]]),
                       se = c(all$se[[3]],
                              rural$se[[3]],urban$se[[3]]),
                       ci_u = c(all$ci[[3,1]],
                                rural$ci[[3,1]],
                                urban$ci[[3,1]]),
                       ci_l = c(all$ci[[3,2]],
                                rural$ci[[3,2]],
                                urban$ci[[3,2]]))

incumbency = ggplot(late.plot,aes(group=group,x=group)) +
  geom_point(aes(y=late*100,x=group),size=2) +
  geom_hline(yintercept = 0,linetype = 'dashed') +
  theme_minimal() +
  ylab('Probability Win Next Election (in p.p.)') +
  xlab('Group of Municipalities') +
  geom_errorbar(aes(ymax = ci_u*100, ymin = ci_l*100), width=0,size=.5) 

#ggsave('incumbency_controls.pdf', plot = incumbency, device = 'pdf',height = 10, width = 10, units = 'cm')

#Figure in appendix: scaterplot ----

rd <- subset(retro, (year <2016 & year>2000) & ineligible2 == 0)

rdA = rd[ruralPop >= (median(rd$ruralPop,na.rm=T))] 
rdB = rd[ruralPop < (median(rd$ruralPop,na.rm=T))] 

sequence=seq(0.01, .15, .01)
xlim = c(0,15)
ylim= c(-.4,.2)
xlimS = c(-.15,.15)
ylimS= c(0,.5)

tt=rdA
tt[,dv:= f_elected]

data.plot.1 <- data.frame(coef = as.numeric(),SD = as.numeric(),n = as.numeric(), stringsAsFactors = FALSE)
for(i in sequence){
  rd.data.plot.1 <- tt[abs(tt$margin) <= i, ]
  weights <- kernelwts(rd.data.plot.1$margin,center = 0,bw = i, kernel="triangular")
  data.plot.1reg <- with(rd.data.plot.1, lm(I(dv) ~ elected*margin+as.factor(year)), weights = weights)
  data.plot.1reg.SE <- sqrt(diag(cluster.vcov(data.plot.1reg, cbind(rd.data.plot.1$muni_code))))
  tempo.a <- data.frame(coef = data.plot.1reg$coefficients[2], SD = data.plot.1reg.SE[2], n = data.plot.1reg$df)
  data.plot.1 <- rbind(data.plot.1,tempo.a)   
}


data.plot.1$cuts <- sequence
data.plot.1$confidence <- with(data.plot.1, qt(.975, df = n - 1) * SD)
data.plot.1$upper <- with(data.plot.1, coef + confidence)
data.plot.1$lower <- with(data.plot.1, coef - confidence)

data.plot.1$description <- 'Probability Winning'
data.plot <- rbind(data.plot.1)

plot.above <- ggplot(data.plot, aes(x = cuts * 100, y = coef, group = description)) +
  geom_line(aes(color = description), size = .5, alpha = .7) +
  geom_line(mapping = aes(y = upper, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_line(mapping = aes(y = lower, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_point(data = data.plot, mapping= aes(x=cuts * 100, y=coef),alpha=.9 ,na.rm=T,size=2) +
  geom_hline(yintercept=0, linetype="dotted",colour='gray10') +
  geom_vline(xintercept=7.91,colour='gray10', alpha=.5) +
  theme_minimal() +
  xlab("") +
  ylab('(2.B)\nLATE P(win, t+1)') +
  scale_color_grey(name="") +
  coord_cartesian(ylim=ylim,xlim=xlim) +
  theme(legend.position="none") +
  geom_text(aes(y = 0.1,label = n),vjust=-.50,size=2,hjust=-.1,col='gray40',angle=60)


dos <- subset(tt, margin > 0)
uno <- subset(tt, margin < 0)
tres <- rbind(uno,dos)

y = (tres$dv)

cut <- cut(tres$margin,400, include.lowest = TRUE)
tmp <- aggregate(y, by=list(cut = cut), FUN=mean, na.rm=T)
tmp1 <- aggregate((tres$margin), by=list(cut = cut), FUN=mean, na.rm=T)
data <- data.frame(margin = tmp1$x, y = tmp$x)

scatterAbove = ggplot() +
  stat_smooth(data = tres, aes(margin,y,group=elected), method="loess",na.rm=T,se=FALSE, level = .95, fill = 'gray80',colour="black", size=1.25) +
  geom_point(data = data[data$margin <0,], aes(margin, y),alpha=.9 ,na.rm=T,size=2) +
  geom_point(data = data[data$margin >0,], aes(margin, y),alpha=.9 ,na.rm=T,size=2,shape = 18) +
  theme_minimal() +
  theme(plot.margin=unit(c(1,0,0,1),"mm")) +
  xlab("") +
  geom_vline(xintercept=0, linetype="dotted") +
  coord_cartesian(ylim=ylimS,xlim=xlimS) +
  scale_y_continuous(labels=scales::percent) +
  scale_x_continuous(labels=scales::percent) +
  ylab("(1.A)\nProbability Win Next Election")

#pdf size: 14 x 7. name: rdd_gov.pdf folder: wp/party presence/latex
allA = arrangeGrob(scatterAbove,plot.above,nrow=1,top = '(2) Rural Municipalities')

tt=rdB
tt[,dv:= f_elected]

data.plot.1 <- data.frame(coef = as.numeric(),SD = as.numeric(),n = as.numeric(), stringsAsFactors = FALSE)
for(i in sequence){
  rd.data.plot.1 <- tt[abs(tt$margin) <= i, ]
  weights <- kernelwts(rd.data.plot.1$margin,center = 0,bw = i, kernel="triangular")
  data.plot.1reg <- with(rd.data.plot.1, lm(I(dv) ~ elected*margin+as.factor(year)), weights = weights)
  data.plot.1reg.SE <- sqrt(diag(cluster.vcov(data.plot.1reg, cbind(rd.data.plot.1$muni_code))))
  tempo.a <- data.frame(coef = data.plot.1reg$coefficients[2], SD = data.plot.1reg.SE[2], n = data.plot.1reg$df)
  data.plot.1 <- rbind(data.plot.1,tempo.a)   
}

data.plot.1$cuts <- sequence
data.plot.1$confidence <- with(data.plot.1, qt(.975, df = n - 1) * SD)
data.plot.1$upper <- with(data.plot.1, coef + confidence)
data.plot.1$lower <- with(data.plot.1, coef - confidence)


data.plot.1$description <- 'Probability Winning'
data.plot <- rbind(data.plot.1)

plot.below <- ggplot(data.plot, aes(x = cuts * 100, y = coef, group = description)) +
  geom_line(aes(color = description), size = .5, alpha = .7) +
  geom_line(mapping = aes(y = upper, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_line(mapping = aes(y = lower, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_point(data = data.plot, mapping= aes(x=cuts * 100, y=coef),alpha=.9 ,na.rm=T,size=2) +
  geom_hline(yintercept=0, linetype="dotted",colour='gray10') +
  geom_vline(xintercept=7.91, colour='gray10',alpha=.5) +
  theme_minimal() +
  xlab("Absolute Distance\nTreatment and Control (p.p.)") +
  ylab('(2.B)\nLATE P(win, t+1)') +
  scale_color_grey(name="") +
  coord_cartesian(ylim=ylim,xlim=xlim) +
  theme(legend.position="none") +
  geom_text(aes(y = 0.1,label = n),vjust=-.50,size=2,hjust=-.1,col='gray40',angle=60)


dos <- subset(tt, margin > 0)
uno <- subset(tt, margin < 0)
tres <- rbind(uno,dos)

y = (tres$dv)

cut <- cut(tres$margin,400, include.lowest = TRUE)
tmp <- aggregate(y, by=list(cut = cut), FUN=mean, na.rm=T)
tmp1 <- aggregate((tres$margin), by=list(cut = cut), FUN=mean, na.rm=T)
data <- data.frame(margin = tmp1$x, y = tmp$x)

scatterBelow = ggplot() +
  stat_smooth(data = tres, aes(margin,y,group=elected), method="loess",na.rm=T,se=FALSE, level = .95, fill = 'gray80',colour="black", size=1.25) +
  geom_point(data = data[data$margin <0,], aes(margin, y),alpha=.9 ,na.rm=T,size=2) +
  geom_point(data = data[data$margin >0,], aes(margin, y),alpha=.9 ,na.rm=T,size=2,shape = 18) +
  theme_minimal() +
  theme(plot.margin=unit(c(1,0,0,1),"mm")) +
  #geom_hline(yintercept=0, linetype="dotted") +
  xlab("Margin") +
  #ggtitle(expression(atop("Discontinuity", atop(italic(""), "")))) +
  geom_vline(xintercept=0, linetype="dotted") +
  coord_cartesian(ylim=ylimS,xlim=xlimS) +
  scale_y_continuous(labels=scales::percent) +
  scale_x_continuous(labels=scales::percent) +
  ylab("(2.A)\nProbability Win Next Election")

#pdf size: 14 x 7. name: rdd_gov.pdf folder: wp/party presence/latex
allB = arrangeGrob(scatterBelow,plot.below,nrow=1,top = '(3) Urban Municipalities')



tt=rd
tt[,dv:= f_elected]

data.plot.1 <- data.frame(coef = as.numeric(),SD = as.numeric(),n = as.numeric(), stringsAsFactors = FALSE)
for(i in sequence){
  rd.data.plot.1 <- tt[abs(tt$margin) <= i, ]
  data.plot.1reg <- with(rd.data.plot.1, lm(I(dv) ~ elected*margin+as.factor(year)))
  data.plot.1reg.SE <- sqrt(diag(cluster.vcov(data.plot.1reg, cbind(rd.data.plot.1$muni_code))))
  tempo.a <- data.frame(coef = data.plot.1reg$coefficients[2], SD = data.plot.1reg.SE[2], n = data.plot.1reg$df)
  data.plot.1 <- rbind(data.plot.1,tempo.a)   
}

data.plot.1$cuts <- sequence
data.plot.1$confidence <- with(data.plot.1, qt(.975, df = n - 1) * SD)
data.plot.1$upper <- with(data.plot.1, coef + confidence)
data.plot.1$lower <- with(data.plot.1, coef - confidence)

data.plot.1$description <- 'Probability Winning - All Municipalities'
data.plot <- rbind(data.plot.1)

plot.above <- ggplot(data.plot, aes(x = cuts * 100, y = coef, group = description)) +
  geom_line(aes(color = description), size = .5, alpha = .7) +
  geom_line(mapping = aes(y = upper, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_line(mapping = aes(y = lower, color = description), lty = "dotted", size = 1, alpha = .7) +
  geom_point(data = data.plot, mapping= aes(x=cuts * 100, y=coef),alpha=.9 ,na.rm=T,size=2) +
  geom_hline(yintercept=0, linetype="dotted",colour='gray10') +
  geom_vline(xintercept=7.91,colour='gray10', alpha=.5) +
  theme_minimal() +
  xlab("") +
  ylab('(2.B)\nLATE P(win, t+1)') +
  scale_color_grey(name="") +
  #ggtitle(expression(atop("All Municipalities", atop(italic(""), "")))) +
  coord_cartesian(ylim=ylim,xlim=xlim) +
  theme(legend.position="none") +
  geom_text(aes(y = 0.1,label = n),vjust=-.50,size=2,hjust=-.1,col='gray40',angle=60)


dos <- subset(tt, margin > 0)
uno <- subset(tt, margin < 0)
tres <- rbind(uno,dos)

y = (tres$dv)

cut <- cut(tres$margin,400, include.lowest = TRUE)
tmp <- aggregate(y, by=list(cut = cut), FUN=mean, na.rm=T)
tmp1 <- aggregate((tres$margin), by=list(cut = cut), FUN=mean, na.rm=T)
data <- data.frame(margin = tmp1$x, y = tmp$x)

scatterAbove = ggplot() +
  stat_smooth(data = tres, aes(margin,y,group=elected), method="loess",na.rm=T,se=FALSE, level = .95, fill = 'gray80',colour="black", size=1.25) +
  geom_point(data = data[data$margin <0,], aes(margin, y),alpha=.9 ,na.rm=T,size=2) +
  geom_point(data = data[data$margin >0,], aes(margin, y),alpha=.9 ,na.rm=T,size=2,shape = 18) +
  theme_minimal() +
  theme(plot.margin=unit(c(1,0,0,1),"mm")) +
  #geom_hline(yintercept=0, linetype="dotted") +
  xlab("") +
  #ggtitle(expression(atop("Discontinuity", atop(italic(""), "")))) +
  geom_vline(xintercept=0, linetype="dotted") +
  coord_cartesian(ylim=ylimS,xlim=xlimS) +
  scale_y_continuous(labels=scales::percent) +
  scale_x_continuous(labels=scales::percent) +
  ylab("(1.A)\nProbability Win Next Election")

#pdf size: 14 x 7. name: rdd_gov.pdf folder: wp/party presence/latex
all = arrangeGrob(scatterAbove,plot.above,nrow=1,top = '(1) All Municipalities')

plot = grid.arrange(all,allA,allB)

#ggsave('scatterplot_incumbency.pdf', plot = plot, device = 'pdf',height = 11.69, width = 8.27, units = 'in')

# Density in the appendix ----

bliss = retro[elected==1 & !is.na(inflation)]
table(bliss$year)

density = ggplot(bliss, aes(inflation)) +
  geom_density() +
  theme_tufte() +
  ylab('Density') +
  xlab('Price Index') +
  coord_cartesian(xlim = c(-20,40)) +
  geom_vline(xintercept = median(bliss$inflation,na.rm=TRUE),size=.25,linetype='dashed') +
  scale_y_continuous(labels = scales::percent)

#ggsave('density.pdf', plot = density, device = 'pdf',width = 12, height = 10, units = 'cm')


# Density tests ----
dens <- subset(retro, (year <2016 & year>2000) & ineligible2 == 0)

margin <- dens$margin
margin_rural <- dens[ruralPop > median(dens$ruralPop, na.rm = TRUE)]$margin
margin_urban <- dens[ruralPop < median(dens$ruralPop, na.rm = TRUE)]$margin

### Plot and save results
rdd <- rddensity(X = margin, vce = 'jackknife',p=1)
plot.density.1 <- rdplotdensity(rdd, X = margin, CIuniform = TRUE,xlabel = 'margin')
#ggsave('density1.pdf', plot = plot.density.1$Estplot, device = 'pdf',height = 10, width = 12, units = 'cm')

rdd <- rddensity(X = margin_rural, vce = 'jackknife',p=1)
plot.density.rural <- rdplotdensity(rdd, X = margin_rural, CIuniform = TRUE,xlabel = 'margin')
#ggsave('density_rural.pdf', plot = plot.density.rural$Estplot, device = 'pdf',height = 10, width = 12, units = 'cm')

rdd <- rddensity(X = margin_urban, vce = 'jackknife', p=1)
plot.density.urban <- rdplotdensity(rdd, X = margin_urban, CIuniform = TRUE,xlabel = 'margin')
#ggsave('density_urban.pdf', plot = plot.density.urban$Estplot, device = 'pdf',height = 10, width = 12, units = 'cm')


# Histogram ----

datum = retro[year == 2008 & elected == 1] # so we don't repeat municipalities

plotdata = datum[!is.na(ruralPop),list(ruralPop)]
plotdata[ruralPop >= median(plotdata$ruralPop,na.rm=T),type := 'Rural']
plotdata[ruralPop < median(plotdata$ruralPop,na.rm=T),type := 'Urban']

plotdata$cutoff <- as.factor((plotdata$ruralPop > median(plotdata$ruralPop,na.rm=T))*1)

plot =  ggplot(aes(x=ruralPop,fill=as.factor(type)),data=plotdata) +
  geom_histogram(binwidth=0.0505,color='white') +theme_tufte() + 
  scale_fill_manual(name='', breaks=c("Rural","Urban"),  values=c('Rural'= 'grey20','Urban'= 'gray80')) +
  xlab(label = 'Share of Rural Population\n(2008)') + ylab(label = 'Number of Municipalities')

#ggsave('histogram.pdf', plot = plot, device = 'pdf',width = 15, height = 12, units = 'cm')